Model of coarsening and vortex formation in vibrated granular rods 
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Neicu and Kudrolli observed experimentally spontaneous formation of the long-range orientational 
order and large-scale vortices in a system of vibrated macroscopic rods. We propose a phenomeno- 
logical theory of this phenomenon, based on a coupled system of equations for local rods density 
and tilt. The density evolution is described by modified Cahn-Hilliard equation, while the tilt is 
described by the Ginzburg-Landau type equation. Our analysis shows that, in accordance to the 
Cahn-Hilliard dynamics, the islands of the ordered phase appear spontaneously and grow due to 
coarsening. The generic vortex solutions of the Ginzburg-Landau equation for the tilt correspond 
to the vortical motion of the rods around the cores which are located near the centers of the islands. 
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Vibrated granular materials exhibit many interesting 
phenomena, including formation of cellular and localized 
patterns, convection, phase separation etc. Re- 
cently, Neicu and Kudrolli studied the dynamics of a 
layer of long cylindrical grains (rods) subjected to ver- 
tical vibration, and discovered a surprising phenomenon 
of spontaneous formation of islands of vertically aligned 
rods which co-exit with the "sea" of randomly packed al- 
most horizontal rods. Subsequently, small islands merge 
and form large islands which typically exhibit collective 
vortical motion of rods. Near the core of the vortex this 
motion has a form of solid body rotation, while farther 
away from the core, the angular velocity decays. 

While the bistability and first order phase transitions 
leading to phase separation and coarsening are typical 
in granular dynamics (observed, for example, in mono- 
layer excitation vibrated powder heaping [^,D, 
electrostatically-driven granular systems Q ) and usually 
caused by the inelasticity of grains, the emergence of the 
vortical motion within the ordered phase is rather unex- 
pected. Experiment shows that rods within a vortex 
are tilted in the azimuthal direction, and slowly drift in 
the direction of the tilt. Authors |^] suggested that the 
drift occurs due to the confinement of the rod vibration 
by its tilted neighbors. 

In this Letter we introduce a continuous phenomeno- 
logical model of the transition to the ordered vorti- 
cal state based on the modified Cahn-Hilliard equation 
governing the dynamics of local rods density and the 
Ginzburg-Landau type equation for the tilt. Our model 
reproduces qualitatively the observed phase separation, 
coarsening and vortex formation. We derive the solutions 
for the stationary vortices and discuss their stability. 

Model. The motion of rods is described by the momen- 
tum conservation equation in the form 
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Here 



is the horizontal velocity of rods, p 



is the density, p is the hydrodynamic pressure, the tilt 
vector n — [rix, riy) is the projection of the rod on x — y 
plane, and n — |n|. The term an/o(n)p, accounts for the 
driving force from the vibrating bottom exerted on the 
tilted rod. According to experiments 0, the driving force 
is proportional to the tilt of rods for small tilt values, but 
the saturates and eventually decays to zero at n > no- 
Finally, the term describes the momentum dissipation 
due to bottom friction. Eq. (^ must be augmented by 
the mass conservation equation in the from 



dtp + div(vp) = 
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In the following we assume that the friction is strong, and 
neglect the inertia term with respect to the friction 
term (^v. Thus, we can express the velocity in the form 



an/o(n)/?) 
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Now, substituting the velocity (||) into the mass conser- 
vation law Eq. (H), we obtain 



dtp = C Miv(Vp-an/o(n)p) 
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To describe phenomenologically the experimentally ob- 
served phase separation and coarsening we employ the 
Cahn-Hilliard approach (see for review |lO| ) . We assume 
that the pressure p can be obtained from the variation of 
certain "free energy" type functional of the p field 
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We adopt the standard form of the free energy taking into 
account the local dynamics and diffusive-type coupling 
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where I is the characteristic length scale related to the rod 
thickness. To account for the bistability and phase sepa- 
ration, function / should have two minima separated by a 



1 



maximum. For simplicity we choose a quartic polynomial 
form of /, and define df /dp = A[p — po){p* — p){5 — p), 
where < po < (5 < p* are characteristic rod densities 
depending on the driving acceleration jll]] . Substituting 
and (||) into Eq. (|^), after appropriate rescaling we 
obtain the modified Cahn-Hilliard equation 



V^p~p{l-p)(5-p) = B, 



(11) 



5tp = -V2 {y^p-p{l-p){5-p)) 
- adiv (n/o(n)(p + po)) , 
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where we keep the same notations for the rescaled vari- 
ables and parameters. 

To close the description we need to add an equation 
for the evolution of tilt n. For density p < 1 the vertical 
orientation of rods corresponding to n = is unstable, as 
rods spontaneously tilt. We assume that the growth rate 
of the instability depends on the rods packing density, so 
we can write for the local dynamics 



dtn = /i(p)n - |n|^n 
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with fi{p) = flo — o-iP, 0,0,1 > some constants. As p in- 
creases, the instability saturates and the equilibrium tilt 
diminishes. In addition, rods interact with each other, 
which leads to the additional spatial derivative operator 
D[n] in Eq.(^). Since the tilt field is not divergence- 
free, from the general symmetry considerations, in the 
lowest (second) order, the "diffusion" operator acting 
on n, takes the form D[n] = f2{p) (Ci V^n + ^2Vdivn) . 
The coefficients ^i_2 in this expression are analogous to 
the first and second viscosity in ordinary fiuids (see [^). 
Function /2 (p) describes the decrease of the spatial cou- 
pling strength as the rods density decreases. We assume 
that in the gas phase (p 0) the spatial coupling be- 
tween the rods is small and their tilt becomes large and 
uncorrelated. Accordingly, we set /2 = p, if p > and 
/2 = otherwise. Finally, we include the simplest term 
describing coupling between the tilt and the density gra- 
dient /3Vp. Combining all these terms we arrive at 
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+ /2(p) (aV^n + aVdivn) + /3Vp. 
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It is convenient to introduce new complex variable 
■0 = "-x + iny. Then, Eq. (^) assumes the form of the 
generalized Ginzburg-Landau equation — + £,2/2) 



dt^={h{p)-\^\^)^ + (]{d^+tdy)p 

+ ./2(p) (evV + |(5. + *ay)V 
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Phase separation. Eq. (|^) exhibits phase separation 
only in a certain range of initial conditions p and param- 
eter S. The total mass conservation yields J J pdxdy = 
5$, where S is the cavity area, and $ is the average 
density determined by the filling fraction. Stationary so- 
lution to Eq. (|7|) obeys 
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where B = const is determined below. In the 
separation regime spatially-homogeneous Eq. 
three real roots pi < p2 < Pa- The root p2 corresponds 
to the unstable solution and pi.3 to the stable ones. Final 
stage of the phase separation leads to one domain of the 
high-density phase (solid) p = p3 of the area Sh and the 
low-density one (gas) p = pi of the area Si — S — Sh- 
From the mass conservation one obtains 



ShP3 + (S - Sh)pi = S^. 



(12) 



Here we neglect the interfacial contributions assuming 
that Si,Sii 3> 1. In addition, Eq. ([l^ ) must be aug- 
mented by the condition that the free energy densities 
of the both phases are equal, which is expressed by the 
relation (so-called area rule, see e.g. Q) 
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[p(l-p)(<5-p)-i3)]dp = 0. 



(13) 



Eq. ( |13| ) and Eq. ( |ll[ ) fix the value of B, and, corre- 
spondingly the roots pi^s. ^From Eq. ( p^ ) it follows that 
the solution with two phases is possible if pi < $ < pa. 
It defines the phase separation region in Fig. ^ 



0.4 



0.2 



gas 






only 








vortices 






gas+clusters 








/ solid 






/ only 



-0.2 0.0 0.2 0.4 0.6 0.8 1.0 1.2 
<I> 

FIG. 1. Phase diagram in $ — (5 plane. Here "I> is related 
to the filling fraction and 5 to the driving acceleration. 

Let US consider radially symmetric vortex solutions to 
Eq. ([lo|). In this case p is a function of the polar ra- 
dius r, and ip can be expressed as tjj = exp(±i6)w{r), 
where w is a complex function, and 9 is the polar an- 
gle (for definiteness we take sign -I-). Using d, 
exp{i6){dr + i/rde), we obtain from Eq. (|l(i| ) 



+ idy = 



dtW = /2(eV^w + ^Vlw*) + fiw - \w\^w + (3pr (14) 



^dr — r ^ is the radial Laplacian 
0, Eq. ([l^ ) possesses a stationary 



where = 9^ + 
operator. For^,/? 
solution in the form w = W{r) exp(i(f)Q) with the real 
positive magnitude W and an arbitrary constant phase 
(j)Q. The terms oc a,/3 still permit the constant phase 
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solutions, but they destroy the continuous phase degen- 
eracy. Indeed, Eq. for /3 = yields 

f2i^+^cos2cl}oMW + hW-W'' = (15) 



and sin2(/)o — 0. Solutions exist only for 0o = 0, tt or 
±7r/2. Function W describes the standard (and well- 
documented) vortex solution to the Ginzburg-Landau 

1/2 

equation with the property W fi for r — + oo 
and W r for r 1 (for its rational approximation 
see e.g. [^). Solutions with = 0,7r describe sinks 
(sources) with zero circulations, whereas the solutions 
with (f>o ~ ±7r/2 are vortices with non-zero circulations. 
The sign of 0o determines the direction of rotation. Near 
the vortex core W ^ r, which corresponds to the solid 
body rotation, as velocity vg ~ W. Far away from the 
core the vortex exhibits differential rotation. 

For Pdrp 0, Eq. (^ has constant phase solutions 
only with (f)Q = 0, tt, i.e. with zero circulation. However, 
it does not guarantee the selection of this solution in the 
bulk of large islands where the density gradient is small. 
It is easy to show that solutions with i^o = are ener- 
getically unfavorable with respect to rotating solutions 
with 00 = ±7r/2 if pdrP is small. Indeed, Eq. (|lj) for 
/3drP = can be written in the form 

^- = -1^ (16) 
with the free energy functional 

dr[f2a\drw\^+r-^\w\^) 

+ ^ {iidr + r-')w*f + CO.) - f,\wf + \w\'] (17) 

Substituting vortex solution w = W{r) exp{i(f>Q) in Eq. 
(|l7|), one obtains after integration (since the calculation 
of the vortex energy is rather straightforward, we refer 
interested readers to Ref. [|l^, p. 11). 

C/ - /i/2(e + 6/2 COS 200 ) log R/ra + const (18) 

where ro ~ 0(1) has the meaning of core radius, and R 
is the outer cutoff radius of integration. As one sees from 
Eq. (|l8|), for physically realizable case 6 > 0, the vor- 
tices with 00 = ±7r/2 have lower energy, and therefore 
are more energetically favorable, and are selected in dy- 
namics. In a general case, the vortex solution would have 
a radius-dependent phase 0o- Near the center where the 
density is almost constant, the phase would be close to 
±7r/2, and near the island border where the density de- 
creases rapidly, the phase should approach or vr. This 
scenario suggests that the azimuthal velocity should grow 
with radius near the core, and decrease near edge of the 
vortex, which is confirmed by our numerical simulations 
(see below) and agrees with experiments. In addition. 



U 



one can show that for /3 > the term /3Vp considered as 
a small perturbation, leads to the drift of the vortex core 
towards the gradient of density, and therefore it stabilizes 
the vortex core near the center of the island. 

Numerical simulations of the Cahn-Hilliard equation 
(^ were performed using an EFT split-step method, and 
the Ginzburg-Landau equation ( p^ was solved using ex- 
plicit method. The domain of integration was 100 x 100 
dimensionless units with periodic boundary conditions, 
number of mesh points/EET harmonics was 256 x 256. 
As initial conditions we used p ^ d with with small ampli- 
tude noise and random initial conditions for tp. Selected 
results are presented in Fig. 2-5. 
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FIG. 2. Evolution of the field \ip\, vortices are shown 
as black dots, white area corresponds to larger \tp\ (and 
p), dark vise versa. Parameters: ^ = 1,^2 = 1, 
(3 = 0.04,/o = 1.5 - W'ji = 1.4 - 0.7p, 
S = 0.3,po ~ 0.25, a = 0.03. Images are shown for 
t = 80(a), 400(6), 1600(c) and t = 3200(d). 

At the initial stage of the evolution many vortices and 
small dense clusters (islands) are created throughout the 
domain of integration (Fig. ^). Islands are seen as 
darker areas on the figure because an increase in density 
p results in the decrease of the amplitude of ip. Some 
islands trap vortices and are practically immobile, others 
don't contain vortices and drift in the direction defined 
by the mean value of the tilt n. With time, small is- 
lands disappear and bigger islands grow (Fig. ||b,c). It 
is interesting to note that due to tilt-driven drift coars- 
ening occurs much faster than in ordinary Cahn-Hilliard 
dynamics. Finally, one big island with the vortex in the 
center is formed (Fig. |^). Corresponding density p and 
the phase field arg?/' are shown in Fig. |^. Surprisingly, 
even far from the island there some "dormant vortices" 
in the low-density phase (gas). These vortices are seen 
as the end-points of the phase singularity lines (lines be- 
tween dark and white) in Fig. Hb. These vortices do not 
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annihilate because in the low-density phase the diffusion 
terms in Eq. (10) are absent. They are not seen in Fig. 
^ because their core size is small in the gas phase. How- 
ever, these vortices can be advected into the high-density 
phase from the edges. 
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FIG. 3. Density p (a), phase arg-i/" 


'y^) and tilt ampli- 



tude \tp\ (c) for t = 4960, other parameters as in Fig. ^ 
Black corresponds to p = = arg-^Z) = — 0, white to 
p = 1.4, argt/j = 2tv, = 1.2. 
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FIG. 4. Velocity field corresponding to Fig. ^. 

The velocity field is calculated using Eq. (||). Fig. ^ 
shows the velocity field corresponding to the calculated 
tilt structure of Fig.^. The rods perform circular motion 
around the vortex center. The azimuthal velocity vs 
r is shown in Fig. ^ The velocity is maximal somewhere 
between the core and the island edge. It qualitatively 
resembles the experimental one 0. Outside the island 
the tilt becomes large and the velocity becomes small. 

In conclusion, we developed a phenomenological model 
of the formation of the vortical ordered state in the sys- 
tem of vertically vibrated rods. Our continuum model 
is based on a Cahn-Hilliard equation for the rods aerial 
density coupled to the Ginzburg-Landau equation for the 
rod tilt. The model reproduces the qualitative features 
of the vortex formation process observed in the recent 
experiment 0|. We thank Toni Neicu and Arshad Ku- 
droUi for stimulating discussions. This work was sup- 
ported by the U.S. Department of Energy under grants 
W-31-109-ENG-38 and DE-FG03-95ER14516. Simula- 
tions were performed at the National Energy Research 
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FIG. 5. Azimuthal velocity ve (solid line) and density p 
(dashed line) vs r for the parameter of Fig. 
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